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TECHNICAL MEMORANDUM , 

UNIFORM RANDOM NUMBER GENERATORS 
INTRODUCTION 


With the increased use of simulation and Monte Carlo methods in all the 
various disciplines of engineering and science, the need for sequences of 
numbers that appear to be drawn from particular probability distributions has 
become increasingly more important. The production of these sequences of 
numbers, or so-called "pseudo- random" numbers, must be simple, fast, 
accurate, and, most importantly, reproducible. The purpose of this paper is 
to present automated procedures for the production of these pseudo- random 
numbers consistent with the above criteria. In automating the random number 
generators, emphasis was placed upon the uniform distribution. 

Most computer library subprogram generators are coded in complex - 
machine language and therefore tend to increase, rather than alleviate, any 
confusion existing on the generation of random variables. Presented herein 
are Fortran subprograms for the generation of uniform pseudo- random numbers 
on the unit interval [0, 1). The generators are of the mixed-multiplicative 
congruential type, and the method is that of Marsaglia and Bray [ 1] . The 
subprograms described are for the Univac 1108, the CDC 3200, and the SDS 
930 digital computers. A method for generating normal random variables from 
uniform variables is also included. 


THEORETICAL CONSIDERATIONS 


There are many methods for generating uniform random variables, 
each having inherent advantages and disadvantages depending upon its utilization 
[2] . The mixed- multiplicative congruential generator is discussed herein. 
Because of varied opinions concerning the appropriateness of this type of 
generator for certain Monte Carlo applications [3-5] , the validity of these 
particular applications will not be addressed. This discussion will consist of 
the basic theory of this type of generator and how it may be used conveniently 
via Fortran subprograms by the method of Marsaglia and Bray [ 1] . 



The congruential method is basically an arithmetric recurrence relation 
involving integers in which each new number is generated from the previous 
number by some deterministic approach. The recurrence relation is initiated 
by some initial value and, at some subsequent point, will redevelop forming a 
closed loop. The length of this closed loop is the period of the generator and, 
hopefully, is nearly equal to the total integer population of the machine, which 
is designated by m. 

The deterministic approach of the congruential generators is the relation, 


X. „ = aX. + c mod (m) (0 ^ X. < m) , 
l+l 1 ' l 


which means that the expression aX^ + c is to be divided by m and X ■ is set 

equal to the remainder. To illustrate this, let m (modulus) = 25, a 
(multiplier) = 7, and c = 1, and let X Q = 3 be the initial value for 


X 1 = 7x3 + 1 mod (25) 

X x = 22 , 

X 2 = 7 x 22 + 1 mod (25) 

X 2 = 5 , 

X = 7x5 + 1 mod (25) 
o 

X 3 = 11 , etc. 


Numbers on [ 0, 1) can be obtained by dividing by m. The method usually is 
called multiplicative if c = 0 and mixed if c * 0. The modulus m is 

normally taken as 2 n for an n-bit binary machine and I0 n for an n- digit 
decimal machine. The constants a and c are chosen to provide speed, a 
long period, and good statistical results [6]. 

A basic problem inherent in congruence method generators is that 
choice of X^, a, and c which will insure a maximum period. The following . 

theorem presented by Hull and Dorbell [2] solves this problem. 
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Theorem: The sequence defined by the congruence relation 


X = aXj + c mod (m) 


has full period m, provided 

1. c is relatively prime to m 

2. a = 1 mod (p) if p is a prime factor of m, 

3. a = mod (4) if 4 is a factor of m. 

Thus, if m is a power of 2, as on a binary machine, we need only to have c 
an odd number and a = 1 mod (4). The proof of this theorem is included in 
Reference [2]. 

The most favorable aspect of the congruence generators is the 
characteristic that a sequence of random digits may be reproduced by simply 
starting the generator with the same initial value. In the paper by Stockmal 
[7] , algorithms are presented that evaluate 

X. = f(i) 
or 


i 


f" 1 (X.) , 


where X. is the ith element of a congruential random number sequence. 

These algorithms can be very useful when only certain parts of a Monte Garlo 
simulation need to be repeated and the generation of the entire sequence of 
random variables is not required. 

Congruence generators have been widely used, and the results have 
been favorable. An extensive list of references may be found in Reference 2. 
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FORTRAN CONGRUENT GENERATORS 


Marsaglia and Bray [ 1] developed a method of generating uniform 
random variables by the congruence method utilizing a single Fortran 
statement. They also describe Fortran programs that mix several such 
generators. 

Initially, consider the handling of integers by the SDS 930 digital 
computer. Here, Fortran integers are stored in 24 bits, and the multiplication 

24 

of 2 integers produces a 24-bit integer mod, 2 . When used in algebraic 

expressions, the sign on an integer I is determined by the relation 


23 


m(I) = 


if 0 ^ . I < 2 


-2 24 + I if 2 23 < I ^ 2 24 - 1 , 


23 23 

which has a range from -2 to 2 - 1. ' We may therefore use the single 

* 23 v ; * 23 ■ 

Fortran statement I = I*A for each random integer on -2 < I < 2 - 1 to 

produce a new random integer thus giving the congruence relation 
24 24 

X = a X^ mod (2 ). Finally dividing by 2 and adding 1/2 will produce 

a uniform variate on [0, 1) , 


U = 1/2 + 1/2 24 . 


The process is moire complicated in the Univac 1108 digital computer 
where l’.s complement arithmetic is. used. .Here, multiplication of two 36- 

? ' * 30 

bit integers' yiei (Is the product mod (2 ) but the sign is handled by ' 


m(I) = 


if 0.5 I < 2 35 


-2 36 + 1 + 1 if I ^ 2 35 
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Therefore, steps must be taken to subtract the integer 1 at certain times to 

36 

represent the proper remainders of 2 . The following instruction, which is a 

correction by Grosenbaugh [ 8] to Marsaglia and Bray’ s original Fortran 
instruction, handles this requirement satisfactorily: 

I = I * K + MINO ( 0, ISIGN (K-l, L) ) 


The two functions MINO and ISIGN are standard Fortran library routines. 

MINO determines the minimum value of a series of integer quantities, and 
ISIGN transfers the sign of the second argument to the absolute value of the 
first argument. The constant K can be chosen for maximum period as shown 
by Van Gel der [6], 

The CDC 3200 digital computer handles integers exactly as does the 
Univac 1108. However, the CDC 3200 stores integers as 24 bits, therefore 
the only change required is the power of 2 in all equations. 

These one-line Fortran generators may be incorporated directly into 
programs easily as they are, but Marsaglia and Bray obtained better statistical 
rc:sults upon combining several generators. As an example, consider the 
following SDS 930 generator: 


L = L * ML 

M = M * MM 

J = 1 + LABS (L)/2 16 

24 

U = 1/2 + (N(J) + L + M]/2 

K = K * MK 

N(J) = K . 



The array N is a 128-element array filled with random numbers previously 
assigned by a one-line generator. In the procedure, J is used to choose from 
the N array; J comes from the random integer L after division by the 
appropriate power of 2. The desired random variable U is formed from the 
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sum of the randomly chosen N array element, the random integer L used to 
find J, and a third additional random integer M. The used element of N is 
replaced by the random integer K. Similar procedures are used for the Univac 
1108 and CDC 3200 generators. 

Subprogram listings are given for the Univac 1108, CDC 3200, and SDS 
930 generators in Appendices A, B, and C, respectively. All three 
subprograms are called in the exact same manner, so that a computer program 
may be converted to run on several machines by simply inserting the 
appropriate generator. Initially, the subprograms must be given a non-zero 
integer to establish the internal 128-cell array. An odd number in the millions 
gives excellent statistical results. The initial call to a generator does not 
produce a useful result. Uniform random variables on the unit interval f 0, 1) 
may be obtained by calling the generator with a zero integer argument. To 
repeat a sequence of variables, the generator is simply reinitialized with the 
same initial value. 


GENERATION OF NORMAL VARIABLES 

As in the case of the uniform pseudo- random number generator, there 
are a variety' of ways to generate normally distributed random variables. A 
method of normal variable generation by Box and Muller [9, 10] is discussed 
here. Their method is simple, fast, and requires very little memory storage. 

The following method is used to generate a pair of random deviates 

(X , X ) from the same normal distribution starting from a pair of uniformly 
x z 

distributed random variables (U , U ) distributed [0, 1). 

1 <u 


X. = (-2 in U ) 2 cos (2 jt U.) 
XX Z 

1 

X = (-2 in U ) 2 sin (2 ir U_) . 

Z i z 


The pair (X , X ) will be normally distributed [ 0, l] and are very reliable in 
X z 

the tails of the distribution. The estimated speed is 6. 01 milliseconds per 
variable on a CDC 3200 digital computer. 
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CONCLUSION 


Methods have been shown and subprograms have been developed for the 
fast, simple, and reproducible generation of pseudo-random numbers for 
uniform and normal probability distributions. The accuracy of the numbers 
must not be assumed to be 100 percent for all utilizations, because potential 
error sources do exist. Marsaglia [5] presents results that indicate that 
every multiplicative generator has a defect that makes it unsuitable for certain 
Monte Carlo applications. Other remote error possibilities are stated in the 
references. However, for most applications the generators presented herein 
give accurate results. 
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APPENDIX A 

FORTRAN LISTING OF UN I VAC 1108 UNIFORM GENERATOR 


FUNCTION RAN(JJ) 



RAN 

10 

. 



RAN 

20 

UNIFORM random number generator for the univac 

1108 ' 


RAN 

30 




RAN 

40 

FIRST CALL MUST BE OF THE FORM X ■ RAN ( J ) » 

WHERE U IS 

AN 

RAN 

50 

INITIAL INTEGER value. subsequent calls MUST 

3E OF'THE 

FORM 

RAN 

60 

X a RAN(O) . 



RAN 

70 




RAN 

80 

DIMENSION N ( 128 ) 



RAN 

90 

Ja J J 



RAN 

100 

IF (J.EO.O) GO TO 2 



RAN 

110 

DO 1 I si # 128 



RAN 

120 

J = J*227157*MINO<0-, IS I GN( 227156, J) ) 



RAN 

130 

N( I )aj 



RAN 

140 

La J J 



RAN 

150 

MaJJ 



RAN 

160 

KaJJ 



RAN 

170 

L = L*2 74693*MlNO(0. I S I GN < 274692 , L ) > 



RAN 

180 

MrM*?<t3l33*MlN0(0. IS IGN( 243132, M) ) 



RAN 

190 

I»1ABS(L)/268435456*1 



RAN 

200 

RAN=.5*FL0AT(N< I )*L*M)*,145519152E-10 



RAN 

210 

KaK*2A9l49*MlN0(0, ISIGN(249148.K) ) 



RAN 

220 

N( I )=K 



RAN 

230 

RETURN 



RAN 

240 

END 



RAN 

250- 
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APPENDIX B 

FORTRAN LISTING OF CDC 3200 UNIFORM GENERATOR 


function ranijji 




RAN 

1G 





RAN 

20 

UNIFORM RANDOM NUMBER GENERATOR 

FOR 

THE CDC 3200 


RAN 

30 





RAN 

40 

FIRST CALL MUST BE OF THE FORM 

X 

s RAN( J) . WHERE 1 J IS 

AN 

RAN 

50 

INITIAL integer value. SUBSEQUENT 

CALLS MUST 8E OF ' THE 

FORM 

RAN 

60 

X = RAn<0) . 




RAN 

70 





RAN 

80 

DIMENSION N<128) 




RAN 

90 

J = JJ 




RAN 

100 

IF (J.EQ.O) 3,1 




RAN 

110 

00 2 1=1.128 




RAN 

120 

J=J*227l57+MlN0 (0, I S I GN < 227156 . J ) ) 




RAN 

130 

N< I )=J 




RAN 

140 

K = JJ 




RAN 

150 

Ms J J 




RAN 

160 

L 3 J J 




RAN 

170 

L=L*274693 + MINO(0, I S I GN < 274692, L ) ) 




RAN 

160 

M«M*243133*MIN0(2, I S I GN < 243132 , M ) ) 




RAN 

190 

IoIABS(L)/65636+l 




RAN 

200 

RAN=.5*FLOAT(N(1)+L*M)*,596046447E 

-07 



RAN 

210 

KsK*249l49*MlNQ<0, IS IGN< 249148, K) ) 




RAN 

220 

N{ I )=K 




RAN 

230 

RETURN 




RAN 

240 

END 




RAN 

250 
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APPENDIX C 

FORTRAN LISTING OF SDS 930 UNIFORM GENERATOR 


FUNCTION RAN(JJ) RAN 10 

RAN 20 

UNIFORM random NUMBER GENERATOR FOR THE SDS 930 ran 30 

RAN 40 

FIRST CALL MUST BE OF THE FORM X ■ RAN ( J ) * WHERE -U IS AN RAN 50 

INITIAL INTEGER VALUE. SUBSEQUENT CALLS MUST BE OF'THE FORM RAN 60 

X = RAN(O) . RAN 70 

RAN BO 

DIMENSION N ( 128 ) RAN 90 

JeJJ RAN 100 

IF (J) 1*3,1 RAN 110 

1 DO 2 I >1*128*1 RAN 120 

JbJ*65539 ran 130 

2 .N ( I ) « J RAN 140 

L»JU RAN 150 

M »JJ RAN 160 

K =jJ RAN 170 

3 L«L*43S7 ran 180 

MoM*9l97 RAN 190 

I«l*lABS(L)/65536 RAN 200 

MANa0.5*FLOAT(N(n*L*M)*0,59604644E-07 RAN 210 

K»K*10757 ran 220 

N(I)bK RAN 230 

RETURN RAN 240 

END RAN 250* 
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